myplot3D <- function(z,
                     x,
                     y,
                     zlab = NULL,
                     xlab = NULL,
                     ylab = NULL,
                     zlim = NULL,
                     xlim = NULL,
                     ylim = NULL,
                     cis = TRUE,
                     nlines = 5,
                     heatmap = NULL,
                     hm.cols = NA,
                     hm.border = NA,
                     plane0 = FALSE,
                     plane.aux = FALSE,
                     plane.aux.pos = NULL,
                     plane.aux.col = NULL,
                     main = NULL,
                     cex.main = 1,
                     ntick = NULL,
                     expand = 0.8,
                     theta = 45,
                     phi = 0,
                     r = sqrt(3),
                     d = 1.2) {
  if (!("plot3D" %in% installed.packages())) {
    warning("Package 'plot3D' required. Installing now.")
    install.packages("plot3D")
  }
  library("plot3D")
  if (cis) {
    if (!(is.array(z) & length(dim(z)) == 3L & dim(z)[3] == 3L)) {
      stop(
        "z: Please supply a three-dimensional array where the third dimension
      provides (1) point estimate, (2) lower bound, and (3) upper bound."
      )
    }
  } else {
    if (!((is.array(z) & length(dim(z)) == 2L | is.matrix(z)))) {
      stop(
        "z: Please supply a two-dimensional array or a matrix."
      )
    }
    z <- array(z, dim = c(dim(z), 1L))
  }
  if (!(is.vector(x) & is.numeric(x) & length(x) == dim(z)[1])) {
    stop("x: Please supply a numeric vector of equal length as the first
         dimension of z.")
  }
  if (!(is.vector(y) & is.numeric(y) & length(y) == dim(z)[2])) {
    stop("y: Please supply a numeric vector of equal length as the second
         dimension of z.")
  }
  if (!(is.null(heatmap)) & !(is.matrix(heatmap))) {
    stop("heatmap: Please supply a matrix.")
  }
  
  length.x <- length(x)
  length.y <- length(y)
  nlines <- nlines
  points.x <- x[round(seq(1, length(x), length.out = nlines))]
  points.y <- y[round(seq(1, length(y), length.out = nlines))]
  if (is.null(ntick)) {
    ntick <- nlines
  }
  if (is.null(zlim)) {
    zlim <- range(z)
  }
  if (!(is.null(heatmap))) {
    bottom.space <- abs(diff(zlim)) / 4
    zlim[1] <- zlim[1] - bottom.space
  }
  if (isTRUE(plane0) & prod(zlim) > 0) {
    if (sum(zlim) < 0) {
      zlim[2] <- 0
    } else {
      zlim[1] <- 0
    }
  }
  if (is.null(xlim)) {
    xlim <- range(x)
  }
  if (is.null(ylim)) {
    ylim <- range(y)
  }
  if (any(is.na(hm.cols))) {
    col <- ramp.col(col = c("gray90", "gray40"),
                    n = 100,
                    alpha = .8)
  } else {
    col <- ramp.col(col = hm.cols, n = 100, alpha = .8)
  }
  
  perspbox(
    z = NULL,
    zlim = zlim,
    zlab = zlab,
    x = NULL,
    xlim = xlim,
    xlab = xlab,
    y = NULL,
    ylim = ylim,
    ylab = ylab,
    main = main,
    cex.main = cex.main,
    bty = 'u',
    col.grid = "white",
    col.panel = "gray95",
    col.axis = "gray80",
    lwd.grid = 2,
    lwd.panel = 2,
    ticktype = "detailed",
    ntick = ntick,
    expand = expand,
    theta = theta,
    phi = phi,
    r = r,
    d = d
  )
  
  ## Bottom Heatmap
  if (!(is.null(heatmap))) {
    length.hm.x <- dim(heatmap)[1] + 1
    length.hm.y <- dim(heatmap)[2] + 1
    hm.x <- seq(min(x), max(x), length.out = length.hm.x)
    hm.y <- seq(min(y), max(y), length.out = length.hm.y)
    hist3D (
      z = (heatmap / max(heatmap)) * bottom.space + zlim[1],
      x = hm.x[-length(hm.x)] + abs(hm.x[1] - hm.x[2]) / 2,
      y = hm.y[-length(hm.y)] + abs(hm.y[1] - hm.y[2]) / 2,
      col = col,
      border = hm.border,
      add = TRUE,
      plot = TRUE,
      colkey = FALSE
    )
    # image3D(x = hm.x,
    #         y = hm.y,
    #         z = zlim[1] + bottom.space,
    #         col = adjustcolor("gray95", alpha.f = 0.1),
    #         border = adjustcolor("white", alpha.f = 0.25),
    #         lwd = 2,
    #         add = T)
  }
  
  
  ## Red 0-Plane
  if (isTRUE(plane0)) {
    p0.x <- seq(min(x), max(x), length.out = ntick)
    p0.y <- seq(min(y), max(y), length.out = ntick)
    image3D(
      x = p0.x,
      y = p0.y,
      z = 0,
      col = adjustcolor("rosybrown1", alpha.f = 0.1),
      border = adjustcolor("white", alpha.f = 0.5),
      lwd = 2,
      add = T
    )
  }
  
  ## Horizontal auxiliaray reference plane
  if (isTRUE(plane.aux)) {
    if (is.null(plane.aux.pos)) {
      stop("plane.aux.pos: Please supply a z-value for the vertical position.")
    }
    if (is.null(plane.aux.col)) {
      plane.aux.col <- adjustcolor("rosybrown1", alpha.f = 0.1)
    } else {
      plane.aux.col <- adjustcolor(plane.aux.col, alpha.f = 0.1)
    }
    p.aux.x <- seq(min(x), max(x), length.out = ntick)
    p.aux.y <- seq(min(y), max(y), length.out = ntick)
    image3D(
      x = p.aux.x,
      y = p.aux.y,
      z = plane.aux.pos,
      col = plane.aux.col,
      border = adjustcolor("white", alpha.f = 0.5),
      lwd = 2,
      add = T
    )
  }
  
  
  ## Draw lines along first dimension
  for (k in 1:nlines) {
    xx <- rep(points.x[k], length.x)
    yy <- y
    zz <- z[which(x == points.x[k]), , 1]
    scatter3D(
      z = zz,
      x = xx,
      y = yy,
      col = "gray40",
      lwd = 2,
      lty = 1,
      type = 'l',
      add = T
    )
    if (cis) {
      for (p in 2:3) {
        xx <- rep(min(x), length.x)
        zz <- z[1, , p]
        scatter3D(
          z = zz,
          x = xx,
          y = yy,
          col = adjustcolor("gray80", alpha.f = 0.04),
          lwd = 0.5,
          lty = 1,
          type = 'l',
          add = T
        )
        xx <- rep(max(x), length.x)
        zz <- z[length(x), , p]
        scatter3D(
          z = zz,
          x = xx,
          y = yy,
          col = adjustcolor("gray80", alpha.f = 0.04),
          lwd = 0.5,
          lty = 1,
          type = 'l',
          add = T
        )
      }
    }
  }
  
  ## Draw lines along second dimension
  for (k in 1:nlines) {
    xx <- x
    yy <- rep(points.y[k], length.y)
    zz <- z[, which(y == points.y[k]), 1]
    scatter3D(
      z = zz,
      x = xx,
      y = yy,
      col = "gray40",
      lwd = 2,
      lty = 1,
      type = 'l',
      add = T
    )
    if (cis) {
      for (p in 2:3) {
        yy <- rep(min(y), length.y)
        zz <- z[, 1, p]
        scatter3D(
          z = zz,
          x = xx,
          y = yy,
          col = adjustcolor("gray80", alpha.f = 0.04),
          lwd = 0.5,
          lty = 1,
          type = 'l',
          add = T
        )
        yy <- rep(max(y), length.y)
        zz <- z[, length(y), p]
        scatter3D(
          z = zz,
          x = xx,
          y = yy,
          col = adjustcolor("gray80", alpha.f = 0.04),
          lwd = 0.5,
          lty = 1,
          type = 'l',
          add = T
        )
      }
    }
  }
  if (cis) {
    persp3D(x,
            y,
            z[, , 2],
            col = adjustcolor("gray40", alpha.f = 0.07),
            add = T)
    persp3D(x,
            y,
            z[, , 3],
            col = adjustcolor("gray40", alpha.f = 0.07),
            add = T)
  }
}